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Abstract 


This  technical  note  provides  a  summary  of  work  carried  out  to  improve,  verify  and  validate 
the  Trident  UNDEX  Interface  computer  code  for  DRDC.  The  improvements  included  the 
incorporation  of  a  capability  to  automatically  generate,  with  minimum  input  from  the  user, 
input  files  for  Vast-USA,  and  a  state-of-the-art  survey  of  potential  flow  methods  for 
calculating  the  dynamics  of  near-field  UNDEX  gas  bubbles.  The  verification  and  validation 
focused  on  the  recently  developed  tightly  coupled  Vast-USA  interface.  In  total  six  different 
problems  were  investigated  for:  (1)  comparison  with  the  loosely  coupled  Vast-USA  Interface; 
(2)  comparison  with  LS-DYNA;  (3)  comparison  with  experiment;  and  (4)  consistency  with 
units.  Recommendations  are  made  for:  (1)  further  improvements  to  the  Trident  UNDEX 
Interface,  including  the  incorporation  of  a  boundary  element  technique  for  modeling  near¬ 
field  UNDEX  bubble  dynamics;  and  (2)  further  verification  and  validation  studies. 

Resume 

La  presente  note  technique  resume  les  travaux  qui  ont  ete  menes  afin  d’ameliorer,  de  verifier 
et  de  valider  le  code  informatique  de  1’ interface  Trident  UNDEX  pour  le  compte  de  RDDC. 
Les  ameliorations  visaient  a  integrer  une  fonction  permettant  de  gerer  automatiquement  les 
fichiers  d’entree  pour  Vast-USA  et  exigeant  un  minimum  d’ intervention  de  la  part  de 
l’utilisateur,  ainsi  qu’un  sondage  a  la  fine  pointe  de  la  technologie  concemant  les  methodes 
d’ecoulement  possible,  et  ce,  pour  calculer  la  dynamique  des  bulles  de  gaz  UNDEX  de  champ 
proche.  La  verification  et  la  validation  etaient  axees  sur  l’interface  Vast-USA,  qui  a  ete 
developpee  recemment  et  est  liee  etroitement  a  1’ interface  UNDEX.  Au  total,  six  problemes 
differents  ont  ete  examines,  soit :  (1)  comparaison  avec  l’interface  Vast-USA  qui  est  liee  de 
maniere  lache;  (2)  comparaison  avec  LS-DYNA;  (3)  comparaison  avec  1’ experience; 

(4)  coherence  avec  les  unites.  Les  recommandations  visent  les  points  suivants  : 

(1)  ameliorations  approfondies  de  T  interface  Trident  UNDEX,  y  compris  integration  d’une 
technique  liee  a  des  elements  de  contour  en  vue  de  modeliser  la  dynamique  des  bulles  de  gaz 
UNDEX  de  champ  proche;  (2)  etudes  de  verification  et  de  validation  approfondies. 
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Executive  summary 


Introduction 

Simulation  of  underwater  explosion  effects  on  ships  and  submarines  using  three-dimensional 
computer  modeling  methods  is  an  indispensable  tool  for  naval  platform  survivability.  At 
DRDC  Atlantic,  this  capability  is  provided  in  part  through  the  use  of  the  commercial  software 
Underwater  Shock  Analysis  (USA).  As  USA  does  not  possess  a  graphical  interface  of  its  own, 
special  software  has  been  developed  in  the  past,  with  funding  from  DRDC,  for  accessing  USA 
through  Martec  Limited’s  Trident  finite  element  analysis  (FEA)  software. 

Results 

The  report  describes  recent  improvements  to  the  Trident  interface  to  USA,  including  the 
development  of  a  closely  coupled  link  between  Trident  and  USA.  A  validation  study,  which 
compares  simulations  obtained  with  the  old  and  new  versions  of  the  Trident/USA  interface 
and  by  employing  different  systems  of  units  in  the  modeling  input,  indicates  that  the  closely 
coupled  link  is  the  correct  functioning.  A  literature  review  of  boundary  element  modeling  of 
gas  bubble  dynamics  is  also  provided. 

Significance 

The  closely  coupled  link  between  Trident  and  USA  allows  simulations  of  explosion  effects  to 
be  performed  with  improved  realism  and  efficiency.  More  realistic  simulations  are  possible 
because  structurally  non-linear  behaviour  can  now  be  included  in  ship  and  submarine 
response.  Furthermore,  efficiency  improvements  allow  simulations  to  be  performed  in  less 
time  and  permit  use  of  more  detailed  structural  models.  The  boundary  element  approach  to 
simulating  underwater  explosion  gas  bubbles  is  of  interest  as  it  may  be  considerably  less 
costly  than  other  viable  methods  for  modeling  close-proximity  bubbles. 

Future  Plans 

Simulations  of  explosion  effects  are  an  important  part  of  survivability  assessments  being 
performed  for  the  Maritime  Force  Protection  (MFP)  Technology  Demonstration  Project 
(TDP).  These  assessments  will  use  detailed,  whole-ship  models  and  will  take  advantage  of  the 
new  closely  coupled  Trident/USA  link.  The  report  makes  a  number  of  recommendations  for 
further  enhancing  the  interface  software,  and  for  implementing  a  boundary-element  based 
capability  for  modeling  gas  bubbles. 
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Sommaire 


Introduction 

La  simulation  des  effets  d’ explosions  sous-marines  sur  les  navires  et  les  sous-marins  a  l’aide 
de  methodes  de  modelisation  informatique  tridimensionnelle  est  indispensable  a  la  survie  des 
plates-formes  maritimes.  A  RDDC  Atlantique,  cette  capacite  est  assuree  au  moyen  du  logiciel 
commercial  Underwater  Shock  Analysis  (USA).  Comme  USA  ne  comporte  pas  sa  propre 
interface  graphique,  un  logiciel  special  a  ete  developpe  dans  le  passe,  a  l’aide  de  fonds  de 
RDDC,  et  ce,  afin  de  permettre  d’acceder  a  USA  au  moyen  du  logiciel  d’ analyse  d’ elements 
finis  (AEF)  Trident  de  Martec  Limited. 

Res  u  I  tats 

Le  rapport  decrit  les  ameliorations  qui  ont  ete  apportees  recemment  a  1’ interface  Trident- 
USA,  y  compris  l’elaboration  d’un  lien  lache  entre  ces  deux  logiciels.  Une  etude  de 
validation,  qui  compare  les  simulations  obtenues  des  ancienne  et  nouvelle  versions  de 
l’interface  Trident-USA,  et  l’emploi  de  systemes  differents  d’unites  pour  la  modelisation 
indiquent  que  le  lien  lache  assure  le  mode  de  fonctionnement  approprie.  Le  rapport  examine 
aussi  la  modelisation  des  elements  de  contour  de  la  dynamique  des  bulles  de  gaz. 

Pertinence 

Le  lien  lache  entre  Trident  et  USA  permet  de  simuler  selon  une  efficacite  et  un  realisme 
accrus  les  effets  d’ explosions.  D’autres  simulations  realistes  sont  possibles,  car  le 
comportement  non  lineaire  sur  le  plan  structural  peut  maintenant  etre  ajoute  a  la  reponse  des 
navires  et  des  sous-marins.  En  outre,  les  ameliorations  apportes  sur  le  plan  de  1’ efficacite 
permettre  de  mener  les  simulations  plus  rapidement  et  de  recourir  a  des  modeles  structuraux 
plus  detailles.  La  demarche  faisant  appel  a  des  elements  de  contour  pour  les  simulations 
d’ explosions  sous-marines  de  bulles  de  gaz  est  interessante  dans  la  mesure  ou  elle  peut 
s’averer  nettement  moins  onereuse  que  d’autres  methodes  viables  de  modelisation  des  bulles  a 
proximite  immediate. 

Plans  futurs 


Les  simulations  d’ effets  d’ explosions  constituent  un  aspect  important  des  evaluations  en 
matiere  de  surviabilite  qui  doivent  etre  menees  dans  le  cadre  du  Projet  de  demonstration  de  la 
technologie  (TDP)  de  la  protection  des  forces  maritimes  (MFP).  Ces  evaluations  feront  appel 
a  des  modeles  detailles  de  navires  complets  et  tireront  profit  du  nouveau  lien  lache  Trident- 
USA.  Le  rapport  fait  etat  d’un  certain  nombre  de  recommandations  visant  T  amelioration 
approfondie  du  logiciel  d’ interface  et  la  mise  en  oeuvre  d’une  capacite  fondee  sur  les  elements 
de  contour  en  vue  de  la  modelisation  de  bulles  de  gaz. 
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1.  Introduction 

This  Technical  Note  (TN-06-10)  provides  a  summary  of  the  work  performed  under  PWGSC 
Contract  No.  W7707-5-3133  entitled  “TRIDENT  UNDEX  INTERFACE”. 

The  ability  to  assess  the  loading,  response  and  survivability  of  naval  structures  and  equipment  to 
withstand  the  effects  of  underwater  explosions  (UNDEX)  is  of  fundamental  importance, 
particularly  for  determining  the  threat  damage  potential  and  developing  effective  counter¬ 
measures.  In  previous  contract  work  with  DRDC  Atlantic  and  DRDC  Suffield,  Martec  Limited 
developed,  and  incorporated  into  it’s  Trident  system,  specialized  software  for  calculating 
UNDEX  loading  on  structures  including  and  interface  to  the  commercial  Underwater  Shock 
Analysis  (USA)  code.  Until  recently,  the  USA  interface  involved  USA  running  in  a  standalone 
mode  (a  loosely  coupled  link).  However,  in  a  recent  contract  with  DRDC  Atlantic  (W7707-04- 
2588),  Martec  developed  a  new  version  of  the  Trident  interface  to  USA  that  allows  Trident’s  Vast 
solver  and  the  USA  code  to  run  together  in  a  tightly  coupled  mode.  In  this  mode,  Vast  performs 
the  solution,  using  loads  generated  by  USA.  The  advantage  of  this  new  version  is  that  it  now 
allows  for  nonlinear  response  analyses. 

The  primary  purpose  of  this  contract  was  to  verify  and  validate  the  new  tightly  coupled  Trident- 
USA  interface  using  a  series  of  test  problems.  This  has  involved  comparing  results  from  linear 
analysis  using  both  loosely  coupled  and  tightly  coupled  versions,  and  UNDEX  analysis  involving 
shock  analysis  using  doubly  asymptotic  approximations  (DAA)  and  bubble  pulse  loading  using 
virtual  mass  approximations  (VMA).  The  verification  and  validation  carried  out  in  this  contract  is 
presented  in  Section  2.0 

The  USA  code  itself  consists  of  several  modules,  each  one  executed  separately,  and  requiring 
separate  input  data  sets.  Preparing  the  latter  is  most  often  tedious  and  error  prone.  Under  this 
contract,  a  capability  to  automatically  generate  these  files  -  with  minimum  input  from  the  user  - 
has  been  implemented  in  the  UNDEX  software.  This  capability  is  described  more  fully  in  Section 
3.0. 

Although  the  specialized  UNDEX  software  developed  by  Martec  offers  a  sufficient  capability  for 
dealing  with  far-field  bubble  dynamics,  such  is  not  the  case  for  the  near-  field.  In  recent  years 
potential  flow  methods  for  calculating  the  dynamics  of  near-field  UNDEX  gas  bubbles  have  been 
investigated.  As  part  of  this  current  contract,  a  survey  of  the  work  in  this  area  has  been  carried 
out,  and  recommendations  made  for  future  implementation  of  such  a  capability  into  the  UNDEX 
software.  The  result  of  this  survey  is  presented  in  Section  4.0 
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2.  Verification  and  Validation  of  the  New  USA  Interface 

In  all  six  different  problems  have  been  considered  for  verification  and  validation.  They  were:  (1) 
the  quarter  sphere  USA  demonstration  problem;  (2)  a  full  sphere  with  diaphragm;  (3)  the  quarter 
cylinder  USA  demonstration  problem;  (4)  the  submerged  cylinder  problem  described  in 
[DeRuntz,  1986];  (5)  the  CPF;  and  (6)  the  Suffield  circular  test  plate  problem. 

2.1  Quarter  Sphere  Analysis 

The  finite  element  model  for  the  quarter  sphere  problem  is  shown  in  Figure  2.1.1.  This  analysis 
uses  a  plane  wave  approximation  (PWA)  to  the  loading.  The  magnitude  of  the  loading  is  such  that 
the  sphere  responds  in  a  linear  manner. 

Computed  velocity  time  histories  for  two  nodes,  134  and  235,  without  damping,  using  both  the 
tightly-coupled  (T-C)  and  loosely-coupled  (L-C)  USA  interfaces,  are  shown  in  Figures  2.1.2  and 
2.1.3. 

Computed  velocity  comparisons  with  mass  proportional  damping  (A  Type)  are  shown  in  Figures 
2.1.4  and  2.1.5.  Computed  velocity  comparisons  with  stiffness  proportional  damping  (B  Type) 
are  shown  in  Figures  2.1.6  and  2.1.7. 

Computed  velocity  comparisons  using  English  units  (inches  and  pounds)  and  metric  units 
(millimeters  and  Newtons)  are  shown  in  Figures  2.1.8  and  2.1.9. 

In  the  tightly-coupled  USA  interface,  there  is  a  Vast  module  for  linear  analysis  (DISP3)  and  a 
Vast  module  for  nonlinear  analysis  (DISP3N).  Computed  velocity  comparisons  using  DISP3  and 
DISP3N  modules  are  shown  in  Figures  2.1.10  and  2.1.11. 

Computed  velocity  comparisons  using  Vast  and  LS-Dyna  are  shown  in  Figures  2.1.12  and  2.1.13. 
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Figure  2.1.1:  Trident  Finite  Element  Model  of  Quarter  Sphere 
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Figure  2.1.2:  T-C  and  L-C  Velocity  Comparison  for  Node  134 
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Figure  2.1.3:  T-C  and  L-C  Velocity  Comparison  for  Node  235 


Ftflffnlut  Fite 


X 

Figure  2.1.4:  T-C  and  L-C  Velocity  Comparison  for  Node  134  -  With  A  Type  Damping 
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Figure  2.1.5:  T-C  and  L-C  Velocity  Comparison  for  Node  235  -  With  A  Type  Damping 
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Figure  2.1.6:  T-C  and  L-C  Velocity  Comparison  for  Node  134  -  With  B  Type  Damping 


DRDC  Atlantic  CR  2006-155 


5 


Fidomal  Fite 


X 

Figure  2.1.7:  T-C  and  L-C  Velocity  Comparison  for  Node  235  -  With  B  Type  Damping 
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Figure  2.1.8:  Velocity  Comparison  for  Node  134  -  English  Units  Versus  Metric  Units 
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Figure  2.1.9:  Velocity  Comparison  for  Node  235  -  English  Units  Versus  Metric  Units 
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Figure  2.1.10:  Velocity  Comparison  for  Node  134  -  Linear  Module  Versus  Nonlinear  Module 
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Figure  2.1.11:  Velocity  Comparison  for  Node  235  -  Linear  Module  Versus  Nonlinear  Module 
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Figure  2.1.12:  Velocity  Comparison  for  Node  134  -  Vast  Versus  LS-Dyna 
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Figure  2.1.13:  Velocity  Comparison  for  Node  235  -  Vast  Versus  LS-Dyna 


2.2  Sphere  with  Diaphragm 

The  quarter  sphere  model  used  in  Section  2.1  had  a  one-on-one  node  and  element  count  for 
structural  and  wetted  surface  elements.  The  sphere  with  diaphragm  model  shown  in  Figure  2.2.1 
represents  the  typical  case  where  there  are  more  structural  elements  than  wetted  surface  elements. 

Computed  velocity  time  histories  for  two  nodes,  134  and  235,  using  both  the  tightly-coupled  (T- 
C)  and  loosely-coupled  (L-C)  USA  interfaces,  are  shown  in  Figures  2.2.2  and  2.2.3. 
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Figure  2.2.1:  Trident  Finite  Element  Model  of  Sphere  With  Diaphragm 
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Figure  2.2.2:  T-C  and  L-C  Velocity  Comparison  for  Node  134 
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Figure  2.2.3:  T-C  and  L-C  Velocity  Comparison  for  Node  235 


2.3  Quarter  Cylinder  Analysis 

The  finite  element  model  for  the  quarter  cylinder  problem  is  shown  in  Figure  2.3.1.  This  analysis 
uses  a  doubly  asymptotic  approximation  (DAA1)  to  a  spherical  shock  wave  loading.  The 
magnitude  of  the  loading  is  such  that  the  cylinder  responds  in  a  linear  manner. 

Computed  velocity  time  histories  for  two  nodes,  1  and  117,  using  both  the  tightly-coupled  (T-C) 
and  loosely-coupled  (L-C)  USA  interfaces,  are  shown  in  Figures  2.3.2  and  2.3.3. 

Computed  velocity  comparisons  using  English  units  (inches  and  pounds)  and  metric  units 
(millimeters  and  Newtons)  are  shown  in  Figures  2.3.4  and  2.3.5. 
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Figure  2.3.1:  Trident  Finite  Element  Model  of  Quarter  Cylinder 
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Figure  2.3.2:  T-C  and  L-C  Velocity  Comparison  for  Node  1 
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Figure  2.3.3:  T-C  and  L-C  Velocity  Comparison  for  Node  1 17 
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Figure  2.3.4:  Velocity  Comparison  for  Node  1  -  English  Units  Versus  Metric  Units 
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Figure  2.3.5:  Velocity  Comparison  for  Node  117  -  English  Units  Versus  Metric  Units 


2.4  Submerged  Cylinder  Analysis 

The  finite  element  model  for  the  submerged  cylinder  problem  is  shown  in  Figures  2.4.1  and  2.4.2. 
This  analysis  uses  a  virtual  mass  approximation  (VMA)  to  a  bubble  pulse  loading.  The  magnitude 
of  the  loading  is  such  that  the  cylinder  responds  in  a  linear  manner.  Two  types  of  boundary 
conditions  were  considered  -  unconstrained  and  pinned. 

Computed  displacement  and  velocity  time  histories  for  the  unconstrained  condition,  for  two 
nodes,  329  and  169,  using  both  the  tightly-coupled  (T-C)  and  loosely-coupled  (L-C)  USA 
interfaces,  are  shown  in  Figures  2.4.3  to  2.4.6. 

Computed  displacement  and  velocity  time  histories  for  the  pinned  condition,  for  node  169,  using 
both  the  tightly-coupled  (T-C)  and  loosely-coupled  (L-C)  USA  interfaces  are  shown  in  Figures 
2.4.7  and  2.4.8. 

Computed  displacement  and  velocity  comparisons  using  English  units  (inches  and  pounds)  and 
metric  units  (millimeters  and  Newtons)  are  shown  in  Figures  2.4.9  and  2.4.10. 
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Figure  2.4.1:  Trident  Finite  Element  Model  of  Submerged  Cylinder 


Figure  2.4.2:  Trident  Finite  Element  Model  of  Submerged  Cylinder  -  Side  View 
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Figure  2.4.3:  T-C  and  L-C  Displacement  Comparison  for  Node  329 
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Figure  2.4.4:  T-C  and  L-C  Velocity  Comparison  for  Node  329 
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Figure  2.4.5:  T-C  and  L-C  Displacement  Comparison  for  Node  169 
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Figure  2.4.6:  T-C  and  L-C  Velocity  Comparison  for  Node  169 
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Figure  2.4.7:  T-C  and  L-C  Displacement  Comparison  for  Node  169  -  Pinned 
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Figure  2.4.8:  T-C  and  L-C  Velocity  Comparison  for  Node  169  -  Pinned 
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Figure  2.4.9:  Displacement  Comparison  for  Node  169  -  English  Units  Versus  Metric  Units 
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Figure  2.4.10:  Velocity  Comparison  for  Node  169  -  English  Units  Versus  Metric  Units 


2.5  CPF  Analysis 


The  finite  element  model  for  the  CPF  is  shown  in  Figures  2.5.1  to  2.5.3.  The  CPF  analyses  have 
been  performed  using  both  a  doubly  asymptotic  approximation  (DAA1)  to  a  spherical  shock 
wave  loading  and  a  virtual  mass  approximation  (VMA)  to  a  bubble  pulse  loading.  The  magnitude 
of  the  loadings  is  such  that  the  vessel  responds  in  a  linear  manner.  Two  types  of  boundary 
conditions  were  considered  -  unconstrained  and  pinned. 
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Computed  displacement  and  velocity  time  histories  for  the  unconstrained  condition,  for  node 
3920,  using  both  the  tightly-coupled  (T-C)  and  loosely-coupled  (L-C)  USA  interfaces,  are  shown 
in  Figures  2.5.4  to  2.5.5,  for  the  bubble  loading,  and  Figures  2.5.6  and  2.5.7  for  the  shock 
loading.. 

Computed  displacement  and  velocity  time  histories  for  the  pinned  condition,  for  node  3920,  using 
both  the  tightly-coupled  (T-C)  and  loosely-coupled  (L-C)  USA  interfaces,  are  shown  in  Figures 
2.5.8  to  2.5.9,  for  the  bubble  loading,  and  Figures  2.5.10  and  2.5.11  for  the  shock  loading. 


Figure  2.5.1:  Trident  Finite  Element  Model  of  the  CPF 
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Figure  2.5.2:  Trident  Finite  Element  Model  of  the  CPF  -  End  View 
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Figure  2.5.3:  Trident  Finite  Element  Model  of  the  CPF  -  Side  View 
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Figure  2.5.4:  T-C  and  L-C  Displacement  Comparison  at  Node  3920  -  Bubble  Loading 
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Figure  2.5.5:  T-C  and  L-C  Velocity  Comparison  at  Node  3920  -  Bubble  Loading 
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Figure  2.5.6:  T-C  and  L-C  Displacement  Comparison  at  Node  3920  -  Shock  Loading 
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Figure  2.5.7:  T-C  and  L-C  Velocity  Comparison  at  Node  3920  -  Shock  Loading 
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Figure  2.5.8:  T-C  and  L-C  Displacement  Comparison  at  Node  3920  -  Bubble  Loading  -  Pinned 
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Figure  2.5.9:  T-C  and  L-C  Velocity  Comparison  at  Node  3920  -  Bubble  Loading  -  Pinned 
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Figure  2.5.10:  T-C  and  L-C  Displacement  Comparison  at  Node  3920  -  Shock  Loading  -  Pinned 
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Figure  2.5.1 1:  T-C  and  L-C  Velocity  Comparison  at  Node  3920  -  Shock  Loading  -  Pinned 


2.6  Suffield  Circular  Plate  Analysis 
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The  experimental  setup  for  the  Suffield  circular  plate  problem  is  shown  in  Figure  2.6.1.  The 
analyses  have  been  performed  using  both  a  doubly  asymptotic  approximation  (DAA1)  to  a 
spherical  shock  wave  loading  and  a  virtual  mass  approximation  (VMA)  to  a  bubble  pulse  loading. 
The  magnitudes  of  the  loadings  are  such  that  the  plate  responds  both  linearly  and  nonlinearly. 

The  finite  element  model  for  the  full  model  is  shown  in  Figure  2.6.2.  The  finite  element  model 
making  use  of  quarter  symmetry  is  shown  in  Figure  2.6.3.  The  results  presented  here  are  for  the 
quarter  symmetry  model.  All  USA  analyses  were  performed  using  the  new  tightly-coupled 
interface. 

A  typical  charge  configuration  is  shown  in  Figure  2.6.4.  Two  charge  scenarios  were  considered: 
(1)  Shot  1  -  weight  =  lOOg  C4  (127  g  TNT  equivalent),  standoff  =  7.5  ft.,  depth  =  13.5;  and  (2) 
Shot2  -  weight  =  250  g  C4  (305  g  TNT  equivalent),  standoff  =  4.5  ft.,  depth  =  16.5.  The  water 
surface  was  taken  to  be  21.0  ft.  above  the  target  panel. 

The  displacement  time  history  comparisons  for  USA  versus  experiment  for  the  shot  1  and  2  shock 
loading  cases  for  linear  analysis  with  no  damping  are  shown  in  Figures  2.6.5  and  2.6.6. 

The  displacement  time  history  comparisons  for  USA  versus  experiment  for  the  shot  1  and  2 
bubble  loading  cases  for  linear  analysis  with  no  damping  are  shown  in  Figures  2.6.7  and  2.6.8. 

The  displacement  time  history  comparisons  for  USA  versus  experiment  for  the  shot  1  and  2 
bubble  loading  cases  for  linear  analysis  with  damping  (AD AMP  =  200.0  and  BDAMP  =  0.0)  are 
shown  in  Figures  2.6.9  and  2.6.10. 

The  displacement  time  history  comparisons  for  USA  linear  analysis  versus  nonlinear  analysis  for 
the  shot  1  and  2  bubble  loading  cases  with  no  damping  are  shown  in  Figures  2.6.1 1  and  2.6.12. 

The  displacement  time  history  comparisons  for  USA  nonlinear  analysis  versus  experiment  for  the 
shot  1  and  2  bubble  loading  cases  with  no  damping  are  shown  in  Figures  2.6.13  and  2.6.14.  The 
material  properties  used  in  the  analyses  were  as  follows:  Young’s  Modulus  =  207000.0;  Poisson’s 
Ratio  =  0.3;  Initial  Yield  Stress  =  550.0;  Isotropic  Tangent  Modulus  =  0.0;  and  Kinematic 
Tangent  Modulus  =  2000.0.  Both  geometric  and  material  nonlinearities  were  included. 


26 


DRDC  Atlantic  CR  2006-155 


Improvement  and  Validation  of  the  Trident  UNDEX  Interface 


Martec  Limited  TN-06-10 


Base  Plate  (10  x  1 0  ft) 
Rests  on  Bottom  of  Pond 


Figure  2.6.1:  Circular  Plate  Experimental  Setup 


Figure  2.6.2:  Trident  Finite  Element  Model  of  Circular  Plate 
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Figure  2.6.3:  Trident  Finite  Element  Quarter  Model  of  Circular  Plate 
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Figure  2.6.4:  Typical  Charge  Configuration  for  Circular  Plate 


28 


DRDC  Atlantic  CR  2006-155 


Improvement  and  Validation  of  the  Trident  UNDEX  Interface 


Martec  Limited  TN-06-10 


External  File 


X 

Figure  2.6.5:  USA  Linear  Analysis  Versus  Experiment  -  Center  Plate  Displacement  -  Shot  1 
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Figure  2.6.6:  USA  Linear  Analysis  Versus  Experiment  -  Center  Plate  Displacement  -  Shot  2 
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Figure  2.6.7:  USA  Linear  Analysis  Versus  Experiment  -  Center  Plate  Displacement  -  Shot  1 
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Figure  2.6.8:  USA  Linear  Analysis  Versus  Experiment  -  Center  Plate  Displacement  -  Shot  2 
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Figure  2.6.9:  USA  Linear  Analysis  with  Damping  Versus  Experiment  -  Center  Plate 
Displacement  -  Shot  1  Bubble  Loading 
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Figure  2.6.10:  USA  Linear  Analysis  with  Damping  Versus  Experiment  -  Center  Plate 
Displacement  -  Shot  2  Bubble  Loading 
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Figure  2.6.1 1:  USA  Linear  Analysis  Versus  Nonlinear  Analysis  -  Center  Plate 
Displacement  -  Shot  1  Bubble  Loading 
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Figure  2.6.12:  USA  Linear  Analysis  Versus  Nonlinear  Analysis  -  Center  Plate 
Displacement  -  Shot  2  Bubble  Loading 
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Figure  2.6.13:  USA  Nonlinear  Analysis  Versus  Experiment  -  Center  Plate  Displacement  - 

Shot  1  Bubble  Loading 
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Figure  2.6.14:  USA  Nonlinear  Analysis  Versus  Experiment  -  Center  Plate  Displacement  - 

Shot  2  Bubble  Loading 
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3.  Automatically  Generating  Input  Files  for  USA 
Analysis 

In  the  case  of  loosely-coupled  USA  analysis,  there  are  five  data  sets  required,  one  for  each  of  the 
modules  VASUSA,  FLUMAS,  AUGMAT,  TIMINT  and  USAVAS.  For  closely-coupled  analysis, 
there  are  four  data  sets,  one  for  each  of  the  modules  VASUSA,  FLUMAS,  AUGMAT  and 
TIMINT. 

In  the  Trident-USA  interface,  the  following  file  naming  convention  is  used:  prefixxxxt.imp, 
prefixx ut.out,  prefix  flu.inp,  prefix  flu. out,  prefix  tim.inp,  prefixtim.out,  prefixx ivt.inp  and 
prefix uvt.out.  Here  “prefix”  is  a  user  defined  file  name  prefix.  The  extension  “inp”  refers  to  input 
file  required  by  the  module  and  “out”  refers  to  the  output  file  created  by  the  module. 

The  new  USA  interface  allows  for  the  automatic  generation  of  the  all  four  data  sets,  in  a  single 
step,  when  using  the  tightly-coupled  link.  When  using  the  loosely-coupled  link,  only  the  data  sets 
for  modules  VASUSA  and  USAVAS  can  be  automatically  generated. 

The  files  automatically  generated  for  the  modules  FLUMAS,  AUGMAT  and  TIMINT  are  in  the 
keyword  format  only.  When  providing  the  necessary  input  parameters,  the  user  has  the  option  to 
apply  data  stored  on  a  UNDEX  charge  data  file. 

The  USA  interface  dialogues  for  generating  the  files  prefixxut.inp,  prefixflu.inp,  prefixaug.inp 
and  prefix  tim.inp  for  tightly-coupled  analyses  are  shown  in  Figure  3.1. 

The  USA  interface  dialogues  for  generating  the  prefixx  ut.inp  and  prefixx ivt.inp  for  loosely- 
coupled  analyses  are  shown  in  Figure  3.2. 
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Figure  3.1:  USA  Interface  Dialogues  for  generating  Input  Files  for  Tightly-coupled  Analysis 
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Figure  3.1:  USA  Interface  Dialogues  for  generating  Input  Files  for  Tightly-coupled  Analysis 

(cont’d) 
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Figure  3.1:  USA  Interface  Dialogues  for  generating  Input  Files  for  Tightly-coupled  Analysis 

(cont’d) 
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Figure  3.1:  USA  Interface  Dialogues  for  generating  Input  Files  for  Tightly-coupled  Analysis 

(cont’d) 
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Figure  3.2:  USA  Interface  Dialogues  for  generating  Input  Files  for  Loosely-coupled  Analysis 


38 


DRDC  Atlantic  CR  2006-155 


Improvement  and  Validation  of  the  Trident  UNDEX  Interface 


Martec  Limited  TN-06-10 


4.  Literature  Review  of  BEM  Techniques  for 
Underwater  Explosion  Modeling 

4.1  Introduction 

Underwater  explosion  cavities  have  been  the  subject  of  considerable  study.  Although  reasonable 
numerical  models  have  been  developed  to  approximate  cavities  in  the  far  field  for  infinite  and 
semi-infinite  domains  as  well  as  simple  geometries,  there  has  been  less  success  in  developing  a 
more  general,  near-field  solution.  This  is  partly  due  to  the  complexity  of  modeling  the  involved 
phenomena;  cavities  migrate,  undergo  expansion/contraction  cycles,  merge  and  depending  on  the 
specific  case,  form  jets  that  strike  nearby  boundaries.  In  spite  of  (or  perhaps  because  of)  this 
complexity,  cavitation  dynamics  has  seen  active  research  for  a  number  of  decades  aimed  at 
improving  the  understanding  of  the  phenomena,  since  bubble  formation  and  collapse  has 
important  military,  industrial,  and  biomedical  applications. 

This  section  describes  a  literature  review  into  boundary  element  techniques  for  near-field 
UNDEX  simulation.  Boundary  element  techniques  lie  in  the  middle  of  techniques  for  UNDEX 
simulations  in  terms  of  physical  representation  and  computation  time.  The  rest  of  the 
introduction  discusses  the  three  main  techniques  used  to  model  UNDEX  cavities  in  order  to  put 
boundary  elements  methods  in  context,  while  the  remainder  of  this  review  focuses  on  boundary 
element  methods,  including  their  theoretical  basis,  limitations,  and  implementation  details. 

There  are  essentially  three  separate  categories  of  models  that  have  been  developed  to  handle 
UNDEX  cavities.  The  first  category  is  simplified  models  [Wardlaw  1998]  and  [Moyer  2003]  that 
attempt  to  reproduce  the  effects  of  an  UNDEX  cavity  through  a  highly  simplified  representation. 
These  methods,  although  computationally  efficient,  are  effectively  limited  to  far-field  cases 
where  the  influence  of  the  geometry  of  the  fluid  domain  on  the  bubble  itself  is  negligible. 

The  other  two  categories  have  seen  increased  interest  in  the  past  three  decades  due  to  their 
applicability  to  near-field  cases  where  the  geometry  of  the  boundary  plays  a  significant  role  in  the 
behavior  of  the  cavity.  These  categories  are  Computational  Fluid  Dynamics  (CFD)  and  Boundary 
Element  Methods  (BEMs). 

CFD  approaches  are  the  most  general  of  the  two  and  come  in  varying  degrees  of  complexity,  but 
for  the  simulation  of  UNDEX  cavities  multi-material  Eulerian  solvers  are  most  common  (e.g.  LS- 
DYNA,  AutoDYN,  Chinook).  These  solvers  work  by  discretizing  the  fluid  volume,  and  then 
applying  finite-difference,  finite- volume  or  finite-element  techniques  to  integrate  the  governing 
partial  differential  equations  forward  in  time.  One  of  the  advantages  of  using  an  Eulerian  CFD 
approach  is  that  acoustic  wave  transmission  (e.g.  shock  wave  from  an  underwater  charge  and 
bubble  pulses)  is  included  in  the  model  by  default.  Furthermore,  CFD  approaches  can  be 
arbitrarily  complicated,  including  everything  from  heat/mass  transfer  to  viscosity  and  turbulence. 
However  CFD  approaches  require  significant  computational  resources  and  suffer  from  a  number 
of  numerical  defects  such  as  smearing  the  cavity  interface,  instabilities  due  to  fluid  mixing, 
numerical  diffusion  of  shocks  and  poor  efficiency.  The  poor  efficiency  is  often  due  to  the 
compressible  (wave-transmission)  portions  of  the  simulation  occupying  only  a  very  small  fraction 
of  the  total  runtime,  since  the  acoustic  timescales  are  on  the  order  of  milliseconds  or  less,  while 
the  bubble  timescales  are  on  the  order  of  tenths  of  seconds  to  seconds.  Consequently  a  large 
portion  of  the  simulation  runtime  is  spent  simulating  an  incompressible  (not  acoustic)  problem 
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using  the  limitations  of  acoustic  wave  propagation  schemes.  Furthermore,  the  generality  afforded 
by  CFD  approaches  may  be  overkill  for  many  practical  applications  involving  cavity  dynamics. 

By  comparison,  BEM  approaches  (which  are  the  primary  focus  of  this  review)  simplify  the 
physical  model  of  the  scenario  to  the  point  where  the  governing  equations  can  be  transformed 
from  differential  equations  into  integral  equations.  The  integral  equations  can  be  solved  using 
iterative  linear  system  solvers,  and  yield  (among  other  things)  the  boundary  velocities.  These 
velocities  are  used  to  advance  the  position  of  the  domain  boundary  in  time,  at  which  point  the 
integral  equations  are  evaluated  again,  new  velocities  are  calculated,  and  the  process  repeats  until 
the  simulation  ends.  BEMs  have  the  advantage  of  having  relatively  modest  memory, 
computation,  and  pre/post  processing  requirements,  but  the  drawback  of  requiring  dramatic 
simplifications  to  the  governing  equations  in  order  to  transform  them  into  integral  equations. 

They  cannot,  for  example,  handle  charge  detonation,  acoustic  wave  propagation  or  viscosity, 
among  other  phenomena.  BEMs  also  suffer  from  the  additional  drawback  of  not  scaling  well, 
since  the  system  of  equations  produced  is  dense,  resulting  in  an  0(N3)  solution  time  (for  direct 
solvers),  unless  more  sophisticated  techniques  such  as  multipole  solvers  are  used.  The  suitability 
of  BEM  techniques  for  a  given  application  is  consequently  largely  problem  dependent.  However 
BEM  techniques  can  be  augmented  with  other  methods  in  order  to  provide  a  more  complete 
representation  of  the  problem.  Using  an  Eulerian  method  to  handle  charge  detonation  and  shock 
propagation,  and  a  BEM  to  handle  cavity  growth  and  collapse  is  one  example  where  a  hybrid 
method  might  achieve  the  benefits  of  both  techniques  without  the  drawbacks  of  either. 

This  section  details  the  results  of  a  literature  review  of  Boundary  Element  Methods  as  applied  to 
UNDEX  cavity  simulation. 

4.2  UNDEX  Cavity  BEM  Overview 

Review  of  a  number  of  different  papers  on  the  use  of  BEM  for  UNDEX  cavities  reveals  that  most 
researchers  have  approached  the  problem  in  a  similar  manner.  Where  not  directly  referenced, 
[Blake  1997a]  and  particularly  [Best,  1993]  provide  good  discussion  of  the  issues  involved. 

There  are  several  steps  that  many  of  the  reviewed  papers  mention  that  are  described  in  the 
following  sections. 

4.2.1  Discretization  and  Boundary  Integral  Representation 

The  domain  is  taken  to  be  the  liquid  only.  Assuming  incompressible,  irrotational  and  invisicid 
flow  allows  the  transformation  of  the  governing  equations  into  a  boundary  integral  formulation. 
These  assumptions  relax  the  physical  problem  from  solving  the  Navier-Stokes  equations 
throughout  a  volume  to  solving  the  Laplace  equation  over  the  volume  boundaries. 

In  order  to  solve  the  Laplace  equation  over  the  surface  of  the  object,  potential  flow  is  used. 
Potential  flow  involves  introducing  a  scalar  field  (the  velocity  potential,  hereafter  referred  to  as 
simply  potential),  the  gradient  of  which  is  the  fluid  velocity  at  any  point  in  the  domain.  This 
potential  field  is  discretized  through  a  superposition  solution,  where  a  collection  of  functions 
(singularities)  that  each  satisfy  the  Laplace  equation  are  distributed  over  the  surface  of  the  object 
in  such  a  manner  as  to  enforce  the  fluid  boundary  conditions.  The  sum  of  these  singularities  also 
satisfies  the  Laplace  equation,  and  computational  effort  is  directed  at  finding  the  strength  of  each 
singularity  so  boundary  conditions  are  met  where  required. 

The  actual  discretization  of  the  boundary  is  similar  to  finite  element  techniques  and  is  usually 
performed  by  specifying  the  nodal  values  of  potential  or  normal  derivative  of  potential,  with 
interpolating  elements  to  specify  how  the  nodal  values  combine  at  points  on  the  boundaries  that 
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are  not  at  the  nodes.  All  boundaries  of  the  fluid  are  discreteized;  any  air-water  surfaces  as  well 
all  rigid  surfaces.  In  two  dimensions,  cubic  splines  seem  to  be  the  most  common  element  type, 
while  in  3D,  triangles  are  used.  Some  papers  refer  to  linear  triangles  [Chahine  1998],  while 
others  refer  to  higher  order  elements  [Blake  1997a].  Most  papers  make  reference  to  quadratic  or 
cubic  elements,  however  3DynaFS  (a  commercial  BEM  code  from  DynaFlow)  makes  use  of 
linear  elements  [Chahine]. 

The  solution  is  subject  to  free  surface  and  rigid  (optionally  moving)  boundary  conditions.  The 
zero  normal  flow  boundary  condition  restricts  fluid  velocities  to  have  no  normal  component  (fluid 
cannot  flow  through  solid  boundaries)  while  the  free-surface  boundary  condition  causes  the 
boundaries  to  move  with  the  local  fluid  velocity.  Initial  conditions  are  specified  as  either  the 
normal  derivative  of  the  potential  functions  for  the  rigid  boundaries,  or  as  the  initial  potential  for 
the  free  surfaces.  When  discretized  as  described  above,  a  dense  set  of  linear  equations  is  formed, 
each  equation  of  which  solves  the  unknown  quantity  (either  potential  or  normal  derivative  of 
potential)  at  a  single  node.  These  values  permit  subsequent  calculation  of  fluid  velocities,  which 
are  used  to  advect  the  free-surfaces  in  time. 

Although  the  system  is  dense,  direct  solvers  are  rarely  used  (although  their  0(N3)  solution  time  is 
used  as  an  example  here  for  lack  of  a  representative  iterative  value).  The  most  commonly  used 
solver  seems  to  the  GMRES  solver,  [Grilli],  [Blake  1997a].  The  coefficients  of  the  system  that  is 
formed  is  purely  geometric  and  based  on  the  fluid  boundaries.  Since  it  is  the  unsteady  evolution 
of  the  boundaries  that  is  of  primary  interest,  the  matrices  must  be  reformed  and  solved  at  every 
timestep.  Consequently  there  is  no  advantage  to  fully  inverting  the  system  to  find  subsequent 
solutions  to  different  system  right-hand-sides.  However  the  solutions  themselves  are  continuous 
in  time,  and  because  of  this,  the  solution  from  the  previous  timestep  will  be  a  very  good  starting 
point  to  initialize  the  solver  for  the  current  timestep,  hopefully  allowing  an  iterative  solver  to 
converge  quickly. 

One  issue  arises  with  distributing  and  also  evaluating  the  singularities  on  the  boundary.  The 
singularities  are,  as  the  name  implies,  singular  (either  1/r  or  1/r2)  and  cannot  be  integrated 
analytically  directly.  This  requires  either  mathematical  tricks  (such  as  a  transformation  to  a 
cylindrical  coordinate  system)  or  numerical  quadrature  to  evaluate  the  surface  integrals  required 
to  form  the  linear  system  [Lin,  2006].  Accurate  evaluation  of  these  weakly  singular  integrals  is 
very  important,  since  they  form  the  diagonals  of  the  linear  system  and  as  such  heavily  influence 
the  systems  numerical  conditioning  and  accuracy  of  the  results. 

4.2.2  Bubble  Representation 

The  bubble-contents  are  represented  as  an  isentropic  ideal  gas.  This  allows  computation  of  the 
internal  pressure  of  the  bubble  based  on  the  ratio  of  specific  heats,  the  starting  pressure  and  ratio 
of  current  volume  to  starting  volume.  The  expression  for  pressure  is  substituted  into  the 
Bernoulli  equation  to  obtain  an  equation  for  the  time  derivative  of  potential  at  the  bubble  surface, 
which  serves  to  link  the  pressure  of  the  gas  in  the  bubble  to  the  integral  equations.  This  time 
derivative  is  integrated  in  time  throughout  the  simulation  to  obtain  the  current  potential  at  the 
bubble  surface  for  any  timestep.  The  initial  conditions  are  created  from  a  spherical  bubble  with 
an  initial  radius  and  internal  pressure.  Full  initial  conditions  require  the  potential  to  be  specified 
over  the  surface  of  the  bubble  and  any  other  free-surfaces.  None  of  the  reviewed  papers  describe 
how  this  initial  potential  is  chosen.  One  possibility  may  be  to  work  from  an  initial  sphere  of  the 
same  radius  as  the  bare-charge,  with  an  equal  potential  to  that  of  the  other  rigid  boundaries, 
creating  stagnant  conditions.  The  pressure  inside  the  bubble  would  be  the  pressure  of  the  charge 
following  detonation,  but  prior  to  any  products  expansion. 
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4.2.3  Main  Solver  Loop 

Once  the  domain  is  discretized  and  the  initial  conditions  set,  the  solver  enters  the  main  time  loop. 
This  includes  several  steps  which  are  summarized  below,  and  described  in  more  detail  in 
[Chahine  1998],  [Blake  1997a]  and  [Tong,  ND]. 

1.  The  free-surface  node  positions  and  surface  potentials  are  integrated  in  time  using 
one  of  any  number  of  integration  schemes  (e.g.  Euler  [Chahine  1998],  Runge-Kutta, 
Adams-Bashforth  [Blake  1997a],  etc.). 

2.  The  geometric  integrals  required  to  form  the  boundary  integral  equation  are 
evaluated,  the  equations  assembled,  and  solved  to  yield  new  values  for  the  normal 
derivative  of  potential  on  free-surfaces  and  potential  on  rigid  surfaces. 

3.  From  the  calculated  values  in  step  2,  tangential  derivatives  of  potential  on  free- 
surfaces  are  calculated.  How  this  is  done  varies  slightly  between  the  papers,  and  is 
discussed  later. 

4.  The  normal  and  tangential  derivatives  of  potential  combine  to  give  the  gradient  of 
potential  on  the  free-surface  points,  which  are  the  fluid  velocities  at  those  points.  The 
fluid  velocities  are  used  to  integrate  nodal  positions  on  the  free-surface  and  for 
determining  the  time  derivatives  of  potential  on  the  free  surface. 

5.  The  volume  of  the  bubble  is  calculated  to  determine  the  new  internal  pressure,  this  is 
not  mentioned  in  any  of  the  UNDEX  related  papers,  but  one  algorithm  [Eberly  2003] 
is  discussed  later. 

6.  The  bubble  internal  pressure  and  free-surface  nodal  velocities  are  used  in  Bernoulli’s 
equation  to  determine  the  time  derivative  of  potential  on  the  free-surfaces. 

7.  With  the  time  derivative  of  potential  and  velocities  calculated  on  the  free-surface,  the 
solver  returns  to  step  1 . 

This  loop  continues  either  until  the  simulation  completes,  or  until  a  jet  forms  that  pierces  the 
bubble.  In  the  event  that  the  jet  pierces  the  bubble,  changes  must  be  made  to  the  discretization  by 
either  introducing  additional  surfaces  into  the  discretization  or  by  adding  a  vortex  ring.  This  is 
discussed  in  a  subsequent  section.  Once  the  discretization  is  modified,  the  solver  returns  to  the 
main  loop  and  continues. 

4.3  Inconsistencies  in  the  Literature 

Although  the  general  solver  structure  seems  to  be  the  same  through  all  the  papers,  there  are 
differences  in  the  implementations.  Choices  such  as  element  order,  methods  of  integration  and 
handling  of  jetting  vary  and  are  discussed  in  the  following  sections. 

4.3.1  Tangential  Derivatives 

The  evaluation  of  tangential  derivatives  is  required  to  complete  the  calculation  of  fluid  velocities 
at  the  free  surfaces.  These  are  in  turn  used  to  calculate  time  derivatives  of  potential  and  to 
integrate  node  positions  in  time.  Since  the  two  dimensional  implementations  typically  use  splines 
to  represent  the  bubble  surface,  it  is  straightforward  to  evaluate  tangential  derivatives  analytically 
from  the  splines. 
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In  three-dimensions,  evaluation  of  the  tangential  derivatives  is  somewhat  more  complex.  Several 
of  the  papers  [Best  1993],  [Blake  1997a],  [Tong]  etc.  fit  multiquadric  surfaces  to  the  area 
immediately  surrounding  the  node  at  which  tangential  derivatives  were  being  evaluated. 
Derivative  calculations  are  then  handled  analytically  on  that  surface.  In  a  paper  for  highly  non¬ 
linear  overturning  wave  simulation  [Grilli  2001]  used  cubic  shape  functions  to  allow  analytical 
evaluation  of  tangential  derivatives,  but  in  a  novel  way  that  used  adjacent  element’s  nodes  rather 
than  increase  the  resolution  of  each  element.  This  was  done  to  avoid  having  to  increase  the 
number  of  nodes  in  the  bubble  representation  (due  to  the  NA3  solution  time  issue)  while  still 
being  able  to  analytically  evaluate  tangential  derivatives  at  a  high  order. 

Although  not  mentioned  in  the  reviewed  papers,  tangential  derivatives  could  also  be  evaluated 
directly  from  the  nodal  potential  values  by  projecting  the  nodal  coordinates  onto  a  tangent  plane 
and  applying  the  concepts  analogous  to  the  discrete  differential  operators  presented  in  [Meyer]. 

Another  method  of  handling  would  be  to  use  cubic  elements  that  enforce  Cl  continuity  but 
implemented  in  an  analogous  fashion  to  Catmull-Rom  splines.  These  would  not  introduce  new 
nodes,  but  would  allow  analytic  calculation  of  tangential  gradients. 

4.3.2  Bubble  Jetting 

One  of  the  key  differences  between  the  reviewed  papers  is  the  manner  in  which  jetting  is  handled. 
When  jetting  occurs,  one  side  of  the  bubble  comes  into  contact  with  the  opposite  side. 
Experimentally  this  results  in  a  torus  shaped  bubble,  and  the  same  trends  are  seen  in  the 
numerical  simulations.  However,  at  the  instant  that  the  jet  contacts  the  opposite  side  of  the 
bubble,  the  topology  of  the  domain  changes,  and  this  introduces  issues  related  to  the 
mathematical  model  used  in  the  simulation  and  in  the  discretization  of  the  simulation. 

4.3.2. 1  Mathematical  Issues 

The  mathematical  issue  that  is  created  is  that  the  flowfield  being  represented  transforms  from  one 
in  which  there  is  zero  circulation  to  one  in  which  there  is  a  finite  circulation.  This  occurs  because 
the  two  contacting  surfaces  of  the  bubble  have  different  surface  potentials,  which  causes  a 
potential  jump.  This  potential  jump  causes  a  flowfield  similar  to  the  streamlines  of  a  dipole 
singularity.  This  is  consistent  with  experiment,  but  is  difficult  to  represent  in  simulations. 

There  are  two  primary  ways  of  handling  this  in  the  literature.  The  first  is  to  introduce  a  cut- 
surface  [Best  1993],  [Blake  1997a]  into  the  discretization  across  the  jet  that  contains  a  potential 
jump.  The  magnitude  of  the  potential  jump  is  such  that  it  induces  the  correct  circulation  in  the 
flow,  while  keeping  the  fluid  domain  simply  connected.  This  cut  surface  is  then  advected  using 
the  fluid  velocity. 

The  second  approach  is  to  insert  a  vortex  ring  [Wang  2005]  around  the  jet  (within  the  bubble), 
which  has  a  proscribed  circulation.  This  has  the  same  effect  as  the  above  method.  However  it 
introduces  the  question  of  what  potential  to  apply  at  the  contact  point,  which  is  strictly  speaking, 
double-valued.  One  option  (the  presented  option)  is  to  simply  use  the  mean  value  of  the  two 
potentials  on  either  side.  What  the  physical  implications  of  this  are  unclear  from  the  paper. 

Both  methods  of  approaching  the  mathematical  issues  simplify  the  problem  by  assuming  that  the 
jet  impact  occurs  at  a  point,  generally  the  first  point  of  contact.  The  potential  values  at  this  point 
are  used  to  determine  the  appropriate  jump  condition  or  circulation  values  for  either  the  cut- 
surface  or  the  vortex  ring.  This  is  a  necessary  simplification,  since  in  reality  the  jet  may  contact 
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the  opposing  cavity  wall  at  many  points,  or  over  a  surface.  Either  of  these  scenarios  results  in  a 
multivalued  jump  condition/circulation  value,  which  may  violate  the  flow-irrotationality 
assumption  of  potential  flow  or  complicate  the  remeshing  step  described  in  the  following  section. 

4. 3. 2. 2  Discretization  Issues 

Regardless  of  whether  a  vortex  ring  or  cut-surface  is  used,  there  are  issues  related  to  the 
discretization  of  the  bubble  surface  when  the  jet  contacts  the  opposing  wall.  At  this  time,  the 
mesh  topology  changes  from  spherical  to  toroidal  [Best  1993],  and  some  remeshing  must  occur  in 
order  to  accommodate  this.  In  two  dimensions,  the  two  nodes  involved  in  the  contact  can  be 
removed,  and  the  two  opposing  surfaces  connected.  Most  papers  also  refer  to  a  smoothing  step  in 
the  remeshing  procedure  to  remove  sharp  discontinuities  that  might  introduce  numerical 
problems. 

In  three  dimensions  this  process  becomes  significantly  more  complex.  Unlike  the  axisymmetric 
case,  the  general  3D  case  does  not  guarantee  that  the  contact  point  occurs  at  two  nodes.  This 
combined  with  mesh  variation  makes  the  connection  process  much  more  complex.  None  of  the 
papers  reviewed  described  how  to  accomplish  this,  although  some  made  vague  references  to  fully 
automated  systems  in  three  dimensions  [Hung  2003]. 

Another  discretization  related  issue  is  that  of  mesh  stretching.  The  jet  itself  is  a  relatively  fine 
feature  which  may  be  under-resolved  if  the  mesh  is  coarse.  As  a  result,  the  tip  of  the  jet  will  have 
very  high  curvature,  but  likely  be  represented  by  only  one  or  two  elements,  with  high  angles 
between  them.  This  will  pollute  calculations  of  surface  integrals  and  pose  two  problems,  one  for 
the  accuracy  of  any  solutions,  and  another  for  the  stability  of  the  code  if  the  system  formed  using 
these  inaccurate  values  has  poor  numerical  properties.  Similarly,  splash  films  [Blake  1997a]  and 
other  features  with  high  curvatures  will  have  the  same  problems  albeit  probably  not  to  the  same 
degree  as  the  jet  itself. 

Several  papers  mention  adaptively  refining  axisymmetric  versions  of  the  algorithms  by 
subdividing  any  segment  that  exceeds  a  certain  length.  This  length  would  logically  be  calculated 
to  take  surface  curvature  into  account  to  maintain  accuracy,  but  no  account  of  that  was  mentioned 
in  the  papers  themselves.  In  three-dimensions,  this  problem  will  worsen,  since  more  panels  will 
be  required  to  represent  the  bubble  to  a  given  accuracy.  Assuming  an  N3  solution  time  in  the 
number  of  nodes,  it  is  likely  that  accuracy  requirements  will  be  relaxed  in  any  three-dimensional 
implementation  below  the  standards  for  two-dimensional  implementations,  and  consequently 
problems  related  to  mesh  resolution  would  be  worse.  [Hung  2003]  claims  that  adaptive 
remeshing  can  introduce  instabilities  and  dramatically  affect  simulation  runtime,  and  instead 
prefers  an  Elastic  Mesh  Technique  proposed  by  Wang.  Again,  the  papers  referring  to  automation 
of  the  three-dimensional  implementations  did  not  mention  any  details  of  how  that  remeshing  was 
performed. 

4. 3. 2.3  Conjecture  on  How  to  Implement  3D  Jetting 

The  following  was  not  mentioned  in  any  of  the  papers,  but  present  some  ideas  on  how  3D  bubble 
jetting  might  be  implemented. 

By  assuming  that  the  impact  of  the  jet  occurs  at  a  single  point  on  the  opposing  wall,  the  closest 
node  on  the  opposing  wall  could  be  found.  The  one-ring  of  elements  surrounding  these  nodes 
would  then  be  removed,  leaving  two  holes  of  approximately  equal  vertex  counts.  The  remeshing 
step  would  consist  of  connecting  these  two  rings  together,  possibly  eliminating  one  or  more  extra 
elements  should  the  number  of  vertices  on  each  ring  be  unequal.  The  nodes  that  are  joined  would 
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either  have  double  valued  potentials  or  have  their  potential  replaced  by  the  mean  potential, 
depending  on  which  scheme  from  the  mathematical  issues  section  was  chosen.  The  resulting 
mesh  would  then  be  smoothed  near  areas  of  high  curvature. 

Another  alternative  would  be  to  create  a  volumetric  representation  of  the  bubble  at  the  time  of 
impact,  and  retesselate  using  something  similar  to  the  Marching  Cubes  algorithm.  Mesh 
simplification  and  shape-based  smoothing  would  then  refine  the  mesh  to  have  high  element 
quality  before  proceeding  with  the  simulation. 

4.4  Issues  Unaddressed  by  the  Reviewed  Papers 

This  section  lists  some  of  the  things  that  were  not  mentioned  in  the  papers  reviewed  that  may 
pose  issues  with  implementing  a  BEM  UNDEX  solver. 

4.4.1  Specifying  Initial  Conditions  for  a  Bubble 

The  calculation  of  the  internal  pressure  of  the  bubble  is  necessary  in  order  to  calculate  the  time 
derivative  of  potential  on  the  free-surfaces.  However  none  of  the  papers  mention  a  way  to  do  this 
for  an  arbitrary  closed  polyhedron.  There  are  a  number  of  methods  in  computer  graphics 
literature  on  how  to  accomplish  this,  one  of  which  is  presented  in  [Eberly  2003]. 

4.4.2  Calculation  of  Pressures  at  Boundaries 

Only  [Blake,  1997a]  and  [Chahine  1998]  make  any  reference  to  actually  calculating  fluid 
pressures.  Integrating  the  conservation  of  momentum  equation  yields  Bernoulli’s  equation,  from 
which  a  liquid  pressure  can  be  found.  The  paper  does  not  go  into  any  great  detail  regarding 
where  this  pressure  is  valid,  it  is  presumed  everywhere  in  the  fluid  domain.  Furthermore,  by 
neglecting  viscosity  potential  flow  results  in  zero  drag  for  fully  immersed  shapes.  This 
phenomenon  is  referred  to  as  D’Alembert’s  paradox.  Consequently  the  nature  of  loadings  that 
can  be  expected  is  unclear.  [Blake  1997a]  does  show  experimental  comparisons,  but  none  appear 
to  be  direct  pressure  trace  comparisons. 

4.4.3  Incorporation  of  Cavitation 

None  of  the  reviewed  papers  made  any  reference  to  incorporating  fluid  cavitation  into  the  BEM 
model.  Whether  this  is  a  separate  area  of  research,  trade  secrets,  or  simply  not  done  is  unclear. 
However  cavitation  plays  a  very  important  role  in  the  overall  response  of  structures  to  underwater 
explosions,  and  its  omission  would  almost  certainly  compromise  results,  particularly  in  the 
nearfield  where  structural  motion  would  have  the  effect  of  cutting  off  the  incident  pressure  wave, 
dramatically  reducing  the  impulsive  loads  seen  by  the  structure. 

One  possibility  would  be  to  introduce  a  scattered  wave  type  boundary  condition  at  the  wetted 
surface  of  the  target  that  essentially  acted  as  a  layer  between  the  BEM  boundary  condition  and  the 
finite  element  load.  If  the  calculated  pressure  dropped  below  the  vapor  pressure,  this  boundary 
condition  would  cut  off  at  the  vapor  pressure.  Rigid  boundary  conditions  in  the  BEM  (prescribed 
normal  derivative  of  potential)  could  be  reformulated  easily  to  include  local  structural  motion 
from  the  finite  element  code.  More  sophisticated  methods  might  include  a  DAA  boundary  rather 
than  a  simple  velocity -based  scattered  wave  boundary. 

It  should  be  noted  that  the  more  sophisticated  acoustic  element  type  cavitation  modeling  (e.g. 

CFA  for  the  USA  code)  would  be  difficult  at  best  to  incorporate  into  a  BEM  simulation  due  to  the 
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thickness  required  for  the  acoustic  layer  particularly  if  jetting  and  very  near  field  effects  were 
desired. 

4.5  Conclusions 

Boundary  element  methods  are  quite  capable  of  simulating  underwater  explosion  cavities  to  a 
reasonable  degree  of  accuracy  in  many  situations.  Axisymmetric  implementations  that  take  a 
bubble  from  initial  conditions  through  to  jetting  are  relatively  common,  and  use  one  of  two 
approaches  for  handling  the  mathematical  and  discretization  challenges  that  occur.  The  first 
approach  is  to  introduce  a  cut  across  the  jet  with  a  potential  jump  that  causes  circulation,  and  the 
second  is  to  introduce  a  vortex  ring  that  accomplishes  the  same  thing.  Both  approaches  require 
the  bubble  boundary  to  be  retesselated  and  smoothed  once  the  jet  contacts  the  opposing  side  of 
the  bubble.  The  choice  on  which  approach  to  use  appears  to  be  largely  based  on  who  is  doing  the 
research,  Best  and  followers  advocate  the  domain  cut,  while  Wang  and  his  followers  the  vortex 
ring.  Both  approaches  appear  to  produce  comparable  results,  and  have  the  same  issues  for 
implementation  in  3D,  although  the  vortex  ring  approach  may  simplify  matters  somewhat  by  not 
advecting  an  additional  surface  requiring  refinement  as  the  simulation  progresses. 

By  far  the  most  challenging  aspect  of  implementing  a  jetting  3D  UNDEX  BEM  solver  is  the 
retesselation  required  when  the  jet  contacts  the  opposing  wall.  Determining  which  elements 
should  remain  in  the  mesh,  which  should  be  discarded,  and  how  to  guarantee  element  quality  are 
not  trivial  tasks.  Some  experimentation  may  be  required  to  determine  an  effective  scheme. 

Other  aspects  to  consider  are,  the  choice  of  element  order,  the  method  of  integration  (both 
spatially  and  temporally)  as  well  as  any  coupling  to  an  external  code  and  its  requirements  (e.g. 
cavitation).  Fortunately  there  is  significant  variation  in  the  literature  regarding  these  issues, 
suggesting  that  a  variety  of  approaches  will  work  well.  Similarly,  the  components  of  the  solver 
that  vary  (e.g.  choice  of  element  order,  integration  methods)  tend  to  be  very  segmented,  so 
different  options  may  be  tried  without  impacting  the  overall  structure  of  the  solver. 
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5.  Summary 

The  PWA  loading  analysis  option  in  the  Vast-USA  interface  has  been  tested  using  the  quarter 
sphere  model.  Excellent  agreement  between  the  tightly-coupled  and  loosely-coupled  Vast-USA 
interfaces  has  been  obtained.  This  has  been  shown  to  be  the  case  with  no  damping  and  with 
various  forms  of  damping.  Identical  results  were  obtained  when  using  either  English  or  metric 
units  for  the  model.  Although  the  loading  was  such  that  the  sphere  responds  in  a  linear  manner, 
an  analysis  was  performed  using  the  nonlinear  module  as  well,  the  latter  giving  the  same  results 
as  the  linear  module.  In  addition,  good  agreement  was  obtained  when  comparing  the  Vast-USA 
analysis  with  a  LS-Dyna-USA  analysis  of  the  same  sphere.  The  analysis  of  the  quarter  sphere 
with  the  diaphragm  demonstrated  the  ability  of  the  tightly-coupled  interface  to  handle  cases 
where  there  are  more  structural  elements  than  wetted  surface  elements. 

The  DAA1  loading  analysis  option  in  the  Vast-USA  interface  has  been  tested  using  the  quarter 
cylinder  model.  Very  good  agreement  between  the  tightly-coupled  and  loosely-coupled  Vast- 
USA  interfaces  has  been  obtained.  Identical  results  were  obtained  when  using  either  English  or 
metric  units  for  the  model. 

The  VMA  loading  analysis  option  in  the  Vast-USA  interface  has  been  tested  using  the  submerged 
cylinder  model.  Reasonably  good  agreement  between  the  tightly-coupled  and  loosely-coupled 
Vast-USA  interfaces  has  been  obtained  for  the  case  where  the  model  is  pinned.  The  displacement 
comparison  for  the  unconstrained  case  was  not  as  good,  as  there  appears  to  more  rigid  body 
drifting  in  the  loosely-coupled  results.  Perhaps  some  further  study  is  warranted  here.  Identical 
results  were  obtained  when  using  either  English  or  metric  units  for  the  model. 

The  CPF  analyses  have  demonstrate  both  the  DAA1  and  VMA  loading  analysis  options  in  the 
Vast-USA  interface  for  a  fairly  large  and  complex  model.  Very  good  agreement  between  the 
tightly-coupled  and  loosely-coupled  Vast-USA  interfaces  has  been  obtained.  Here  the  rigid  body 
drifting  in  the  unconstrained  model  is  not  as  pronounced  as  in  the  case  of  the  unconstrained 
submerged  cylinder. 

The  Suffield  circular  plate  analysis  has  demonstrated  the  ability  of  the  Vast-USA  interface  to 
accurately  predict  the  structural  response  of  such  systems  to  both  shock  and  bubble  induced 
pressure  loading.  The  results  show  the  importance  of  including  the  nonlinear  effect,  particularly 
for  the  bubble  pulsating  pressure  loading.  The  effect  of  damping  is  also  evident.  The  predictions 
for  Shot  1  are  closer  to  experimental  results  than  for  Shot  2.  Obviously  there  is  a  lot  more  work 
that  could  be  done  here.  For  example,  damping  could  be  modeled  more  accurately,  which  would 
improve  the  predictions  for  both  shots.  Also,  it  is  possible  that,  in  the  case  of  Shot  2,  the  bubble  is 
too  close  to  the  plate  for  accurate  analysis  using  USA. 

The  USA  interface  dialogues  for  generating  the  input  files  for  tightly-coupled  analyses  has 
proven  to  be  extremely  useful  in  carrying  out  the  various  verification  and  validation  analyses 
described  above.  As  stated,  the  input  files  for  the  USA  modules  are  generated  using  the  keyword 
format.  No  effort  was  made  to  provide  the  same  capability  for  the  loosely-coupled  Vast-USA 
interface.  If  use  of  the  latter  is  to  be  continued,  an  attempt  should  be  made  to  provide  a  similar 
automated  procedure  for  it. 

The  literature  review  of  BEM  techniques  for  UNDEX  modeling  of  bubble  dynamics  has  shown 
that  the  extension  of  the  Trident  UNDEX  software  to  include  such  a  solver  is  both  doable  and 
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desirable.  Since  most  work  to-date  has  centered  on  up-to  jetting  only,  it  is  recommended  that  the 
effort  concentrate  on  this  aspect  of  bubble  dynamics  only.  Following  this,  the  much  more 
challenging  problem  of  implementing  a  jetting  UNDEX  BEM  solver  could  be  tackled. 
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